library (tmaptools)
library (RUMBA)
library (tidyverse)
library (sf)
library (ggmap)
library (osmdata)

1 Geocodificando

Pero… qué es geocodificar?

Es el proceso de asignar coordenadas (ergo, localización espacial) a uno o varios atributos alfanuméricos que constituyen una dirección. Se requiere de 4 argumentos fundamentales para ello: Calle, Altura o número, Localidad y Estado o Provincia. Otros atributos mejoran la calidad y la precisión: Barrio, País, Código Postal, etc.

Para geocodificar hacemos uso de APIs y servicios geográficos, en este caso, veremos cómo usar el de OpenStreetMap o también conocido como OSM (Nominatim)

1.1 OpenStreetMap

OpenStreetMap es un proyecto colaborativo para crear mapas editables y libres. Los mapas se crean utilizando información geográfica capturada con dispositivos GPS móviles, ortofotografías y otras fuentes libres. Esta cartografía, tanto las imágenes creadas como los datos vectoriales almacenados en su base de datos, se distribuye bajo licencia abierta Licencia Abierta de Bases de Datos.

Los contribuidores más entusiastas mapean barrios completos utilizando herramientas GPS para enviar información local completa, actualizada y precisa a OpenStreetMap. Varias empresas y entidades públicas que producen información geográfica también contribuyen al permitir que sus datos sean incluidos. Existen equipos profesionales de contribuidores que se coordinan para agregar y mantener actualizada información georeferenciada de límites políticos, calles, edificios, negocios y otros puntos de interés.

Toda la información disponible en OpenStreetMap puede ser descargada y reutilizada por cualquier persona, ya sea accediendo al mapa online, obteniendo una copia completa de la base de datos, o accediendo a los datos vía API.

Empecemos!

Vamos a trabajar con datos de una base de datos de entros de salud municipales de Cordoba capital.

setwd("C:/Users/Elian/OneDrive/Escritorio/Ciencia de datos/Geo en Salud/Clase 4") #definiendo el directorio de trabajo
Centros_Salud <- readxl::read_xlsx("Listado_Centros_de_Salud_Municipales.xlsx")#deben colocar el directorio en el cual descomprimieron el archivo

#vamos a explorar el objeto
head (Centros_Salud)
## # A tibble: 6 x 8
##       Z    Nº `Centro de Salud`  Calle      Altura Localidad Horarios Telefono  
##   <dbl> <dbl> <chr>              <chr>      <chr>  <chr>     <chr>    <chr>     
## 1     1     1 GENERAL MOSCONI    Pedro Naon 1330   Cordoba   7-21 hs  4335102 /~
## 2     3     2 LOS SAUCES         Algarrobo  351    Cordoba   7-19 hs  4338516   
## 3     2     3 LOS PINOS          Pje. Gabl~ 3786   Cordoba   7-14hs   4339101   
## 4     1     4 NUEVA ITALIA       Jose de Q~ 2430   Cordoba   7-17 hs  4335123   
## 5     2     5 YOFRE  NORTE       Rene de C~ 1875   Cordoba   7-14 hs  4339103   
## 6     6     6 VILLA RIVERA INDA~ Chamico e~ <NA>   Cordoba   7-16.12~ 03543-449~
summary(Centros_Salud)
##        Z              Nº         Centro de Salud       Calle          
##  Min.   :1.00   Min.   :  1.00   Length:100         Length:100        
##  1st Qu.:2.00   1st Qu.: 25.75   Class :character   Class :character  
##  Median :4.00   Median : 51.50   Mode  :character   Mode  :character  
##  Mean   :3.64   Mean   : 51.15                                        
##  3rd Qu.:5.00   3rd Qu.: 76.25                                        
##  Max.   :6.00   Max.   :101.00                                        
##     Altura           Localidad           Horarios           Telefono        
##  Length:100         Length:100         Length:100         Length:100        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
## 
#cuántos valores nulos tenemos en los datos de dirección?
sum(is.na(Centros_Salud$Calle))
## [1] 0
sum(is.na(Centros_Salud$Altura)) #Tenemos varios valores nulos
## [1] 26

Para acceder a las funcionalidades de la API usaremos la libería tmaptools, que tiene una amplia funcionalidad, entre ellas, conectarse a los servicios de OSM.

# primero tenemos que construir el campo de dirección tal como lo requerimos
Centros_Salud_Direccion <- Centros_Salud %>% 
  mutate(DirCompleta = paste(Calle, Altura, ",", Localidad, ", Argentina"))

#vamos a explorar el objeto
head(Centros_Salud_Direccion)
## # A tibble: 6 x 9
##       Z    Nº `Centro de Salu~ Calle Altura Localidad Horarios Telefono
##   <dbl> <dbl> <chr>            <chr> <chr>  <chr>     <chr>    <chr>   
## 1     1     1 GENERAL MOSCONI  Pedr~ 1330   Cordoba   7-21 hs  4335102~
## 2     3     2 LOS SAUCES       Alga~ 351    Cordoba   7-19 hs  4338516 
## 3     2     3 LOS PINOS        Pje.~ 3786   Cordoba   7-14hs   4339101 
## 4     1     4 NUEVA ITALIA     Jose~ 2430   Cordoba   7-17 hs  4335123 
## 5     2     5 YOFRE  NORTE     Rene~ 1875   Cordoba   7-14 hs  4339103 
## 6     6     6 VILLA RIVERA IN~ Cham~ <NA>   Cordoba   7-16.12~ 03543-4~
## # ... with 1 more variable: DirCompleta <chr>
summary(Centros_Salud_Direccion)
##        Z              Nº         Centro de Salud       Calle          
##  Min.   :1.00   Min.   :  1.00   Length:100         Length:100        
##  1st Qu.:2.00   1st Qu.: 25.75   Class :character   Class :character  
##  Median :4.00   Median : 51.50   Mode  :character   Mode  :character  
##  Mean   :3.64   Mean   : 51.15                                        
##  3rd Qu.:5.00   3rd Qu.: 76.25                                        
##  Max.   :6.00   Max.   :101.00                                        
##     Altura           Localidad           Horarios           Telefono        
##  Length:100         Length:100         Length:100         Length:100        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##  DirCompleta       
##  Length:100        
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
#cuántas direcciones tenemos?
nrow (Centros_Salud_Direccion)
## [1] 100

Vamos a geocodificar nuestras direcciones

Centros_Salud_Direccion <- Centros_Salud_Direccion %>% 
  mutate(geo=geocode_OSM(  #me devuelve las coordenadas (geocode_OSM) en un campo adicional a los originales
    Centros_Salud_Direccion$DirCompleta,  
    return.first.only = TRUE,  # me devuelve solo el primer resultado encontrado
    keep.unfound = TRUE,  # que mantenga aquellas direcciones que No pudo geolocalizar, con un null en el campo "geo"
    details = FALSE, # que no devuelva detalles, id de OSM u objetos cercanos
    #as.data.frame = NA,
    as.sf = FALSE, # que lo deje como una tabla y no lo convierta en un spatial frame
    geometry = "point", # que la geometría sea de tipo PUNTO (es decir, un par de coordenadas)
    server = "https://nominatim.openstreetmap.org" # el servicio que deseo usar: Nominatim
  ))

#qué obtuvimos?
head(Centros_Salud_Direccion)
## # A tibble: 6 x 10
##       Z    Nº `Centro de Salu~ Calle Altura Localidad Horarios Telefono
##   <dbl> <dbl> <chr>            <chr> <chr>  <chr>     <chr>    <chr>   
## 1     1     1 GENERAL MOSCONI  Pedr~ 1330   Cordoba   7-21 hs  4335102~
## 2     3     2 LOS SAUCES       Alga~ 351    Cordoba   7-19 hs  4338516 
## 3     2     3 LOS PINOS        Pje.~ 3786   Cordoba   7-14hs   4339101 
## 4     1     4 NUEVA ITALIA     Jose~ 2430   Cordoba   7-17 hs  4335123 
## 5     2     5 YOFRE  NORTE     Rene~ 1875   Cordoba   7-14 hs  4339103 
## 6     6     6 VILLA RIVERA IN~ Cham~ <NA>   Cordoba   7-16.12~ 03543-4~
## # ... with 8 more variables: DirCompleta <chr>, geo$query <chr>, $lat <dbl>,
## #   $lon <dbl>, $lat_min <dbl>, $lat_max <dbl>, $lon_min <dbl>, $lon_max <dbl>
summary (Centros_Salud_Direccion$geo) #vamos a explorar el componente geo del objeto que nos devuelve OSM
##     query                lat              lon            lat_min      
##  Length:100         Min.   :-33.27   Min.   :-64.47   Min.   :-33.27  
##  Class :character   1st Qu.:-31.45   1st Qu.:-64.24   1st Qu.:-31.45  
##  Mode  :character   Median :-31.42   Median :-64.21   Median :-31.42  
##                     Mean   :-31.52   Mean   :-64.13   Mean   :-31.52  
##                     3rd Qu.:-31.38   3rd Qu.:-64.14   3rd Qu.:-31.38  
##                     Max.   :-31.32   Max.   :-62.19   Max.   :-31.32  
##                     NA's   :55       NA's   :55       NA's   :55      
##     lat_max          lon_min          lon_max      
##  Min.   :-33.27   Min.   :-64.47   Min.   :-64.47  
##  1st Qu.:-31.45   1st Qu.:-64.24   1st Qu.:-64.24  
##  Median :-31.42   Median :-64.21   Median :-64.20  
##  Mean   :-31.52   Mean   :-64.13   Mean   :-64.13  
##  3rd Qu.:-31.38   3rd Qu.:-64.15   3rd Qu.:-64.14  
##  Max.   :-31.32   Max.   :-62.19   Max.   :-62.19  
##  NA's   :55       NA's   :55       NA's   :55
# Base final con coordenadas
Centros_Salud_Direccion <- Centros_Salud_Direccion %>% 
  mutate(lon=geo$lon, lat = geo$lat) %>% # el campo geo se arma como lista, tenemos que seleccionar una lon y una lat
  select(`Centro de Salud`, DirCompleta, lon, lat)

#exploremos el objeto final...
summary (Centros_Salud_Direccion)
##  Centro de Salud    DirCompleta             lon              lat        
##  Length:100         Length:100         Min.   :-64.47   Min.   :-33.27  
##  Class :character   Class :character   1st Qu.:-64.24   1st Qu.:-31.45  
##  Mode  :character   Mode  :character   Median :-64.21   Median :-31.42  
##                                        Mean   :-64.13   Mean   :-31.52  
##                                        3rd Qu.:-64.14   3rd Qu.:-31.38  
##                                        Max.   :-62.19   Max.   :-31.32  
##                                        NA's   :55       NA's   :55
#qué nos falta para que sea un objeto geográfico? Agregar el sistema de referencia!
Centros_Salud_Direccion <- Centros_Salud_Direccion %>% 
  filter(!is.na(lon), !is.na(lat)) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)


summary (Centros_Salud_Direccion)
##  Centro de Salud    DirCompleta                 geometry 
##  Length:45          Length:45          POINT        :45  
##  Class :character   Class :character   epsg:4326    : 0  
##  Mode  :character   Mode  :character   +proj=long...: 0
class(Centros_Salud_Direccion)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

1.2 USIG

Hay otros servicios para trabajar datos!

USIG es la Unidad de Sistemas de Información Geográfica del Gobierno de la Ciudad de Buenos Aires y brinda un servicio de geocodificación (entre otras cosas) de datos de la Ciudad de Buenos Aires de manera gratuita!!!

La librería que nos facilita el uso de este servicio se llama RUMBA, la cual la instalamos con el clásico install.packages(“RUMBA”).

RUMBA es un conjunto de herramientas para el análisis de la Región Urbana Metropolitana de Buenos Aires usando R. Incluye funciones que permiten obtener coordenadas precisas (longitud y latitud) que corresponden a direcciones dentro de los límites de la Región Urbana Metropolitana de Buenos Aires o bien de barrios vulnerables de la Ciudad de Buenos Aires.

Las funciones consultan la API del Normalizador de direcciones y de Déficit Habitacional de la USIG. Ademas de las coordenadas, se obtiene la dirección normalizada (escrita de forma inequívoca).

Les proponemos hacer un ejercicio con las postas de vacunación de CABA

Postas_Vacunacion <-  readxl::read_xlsx("postas_vacunacion_covid.xlsx")#deben colocar el directorio en el cual descomprimieron el archivo

#vamos a explorar el objeto
head(Postas_Vacunacion)
## # A tibble: 6 x 7
##   categoria clasif_int   efector     tipo_de_ef       direccion   barrio  comuna
##   <chr>     <chr>        <chr>       <chr>            <chr>       <chr>   <chr> 
## 1 Público   Atención a ~ Oficina Na~ Sistema Público~ Libertad 1~ "RETIR~ Comun~
## 2 Público   Atención a ~ Sede Oscoe~ Sistema Público~ Bernardo d~ "BERNA~ MONSE~
## 3 Público   Atención a ~ Club San L~ Sistema Público~ Av. La Pla~ "BOEDO" Comun~
## 4 Público   Atención a ~ Club Itali~ Sistema Público~ Yerbal 150  "CABAL~ Comun~
## 5 Público   Atención a ~ Parque Roca Sistema Público~ Av. Cnel. ~ "CNEL.~ VILLA~
## 6 Público   Atención a ~ Club Glori~ Sistema Público~ Bragado 68~ "CABA\~ MATAD~
summary(Postas_Vacunacion)
##   categoria          clasif_int          efector           tipo_de_ef       
##  Length:85          Length:85          Length:85          Length:85         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##   direccion            barrio             comuna         
##  Length:85          Length:85          Length:85         
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character
nrow(Postas_Vacunacion)
## [1] 85

El dataframe tiene una columna que hace referencia a la dirección de los centros de vacunación. Pero no tenemos ni latitud ni longitud. Ya sabemos que podemos geocodificar los datos! Vamos a probar el servicio de la USIG.

Con la función mutate_USIG_geocode podemos agregar las columnas de lon y lat y así tener las coordenadas. La dirección debe estar expresada como “calle altura, partido”, “calle altura, municipio”, “calle y calle, partido”, o “calle altura, municipio”. El partido o municipio son opcionales. De no ser aclarados, y encontrarse múltiples direcciones que coincidan con la búsqueda, se entregaran las coordenadas dentro de la Ciudad Autónoma de Buenos Aires (si existieran), o en su defecto las del primer partido -por orden alfabético- donde se haya encontrado la dirección.

En resumen: es mejor incluir partido o municipio en las direcciones a georeferenciar.

Postas_Vacunacion_Direccion <- mutate_USIG_geocode(Postas_Vacunacion, "direccion") #paciencia. Puede tardar un poco! Utilizar la función es tan simple como indicar el DF y el nombre del campo que tiene la dirección que vamos a trabajar.

#vamos a explorar el objeto
head(Postas_Vacunacion_Direccion)
##   categoria                 clasif_int                              efector
## 1   Público Atención a Adultos Mayores               Oficina Nacional Scout
## 2   Público Atención a Adultos Mayores                         Sede Oscoema
## 3   Público Atención a Adultos Mayores Club San Lorenzo (Sede Av. La Plata)
## 4   Público Atención a Adultos Mayores                        Club Italiano
## 5   Público Atención a Adultos Mayores                          Parque Roca
## 6   Público Atención a Adultos Mayores              Club Glorias Argentinas
##                                               tipo_de_ef
## 1 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 2 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 3 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 4 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 5 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 6 Sistema Público de la Ciudad (Posta extrahospitalaria)
##                  direccion           barrio        comuna
## 1            Libertad 1282           RETIRO      Comuna 1
## 2 Bernardo de Irigoyen 330 BERNARDO DE 330"     MONSERRAT
## 3        Av. La Plata 1770            BOEDO      Comuna 5
## 4               Yerbal 150        CABALLITO      Comuna 6
## 5      Av. Cnel. Roca 4200  CNEL. AV. 4200" VILLA SOLDATI
## 6             Bragado 6875            CABA"     MATADEROS
##                address_normalised       lon       lat
## 1             LIBERTAD 1282, CABA -58.38452 -34.59320
## 2 IRIGOYEN, BERNARDO DE 330, CABA -58.38049 -34.61215
## 3         LA PLATA AV. 1770, CABA -58.42459 -34.63538
## 4                YERBAL 150, CABA -58.43292 -34.61538
## 5      ROCA, CNEL. AV. 4200, CABA -58.44747 -34.67549
## 6              BRAGADO 6875, CABA -58.50989 -34.65928
summary(Postas_Vacunacion_Direccion)
##   categoria          clasif_int          efector           tipo_de_ef       
##  Length:85          Length:85          Length:85          Length:85         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##   direccion            barrio             comuna          address_normalised
##  Length:85          Length:85          Length:85          Length:85         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##       lon              lat        
##  Min.   :-58.52   Min.   :-34.68  
##  1st Qu.:-58.47   1st Qu.:-34.63  
##  Median :-58.43   Median :-34.62  
##  Mean   :-58.43   Mean   :-34.62  
##  3rd Qu.:-58.40   3rd Qu.:-34.59  
##  Max.   :-58.36   Max.   :-34.55  
##  NA's   :33       NA's   :33
#Solo nos queda transformarlo a objeto geográfico
Postas_Vacunacion_Direccion <- Postas_Vacunacion_Direccion %>% 
  filter(!is.na(lon), !is.na(lat)) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)


summary (Postas_Vacunacion_Direccion)
##   categoria          clasif_int          efector           tipo_de_ef       
##  Length:52          Length:52          Length:52          Length:52         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##   direccion            barrio             comuna          address_normalised
##  Length:52          Length:52          Length:52          Length:52         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##           geometry 
##  POINT        :52  
##  epsg:4326    : 0  
##  +proj=long...: 0
class(Postas_Vacunacion_Direccion)
## [1] "sf"         "data.frame"

2 Agregando contexto a la información

Pudimos geolocalizar gran parte de nuestras direcciones de interés, aunque ya vimos que no es un proceso infalible… pero cuando trabajamos con datos gegráficos es importante poder tener una referencia, ubicar los puntos en contexto, pero ¿Qué pasa cuando no tenemos un mapa?

Tenemos una opción! Podemos agregar capas base, que sería como agregar “una foto” sobre la cual posicionaremos los datos de interés.

El primer paso para obtener las capas bases que necesitamos es obtener el bounding box del mapa que queremos descargar. ¿Qué utilizamos?… así es! getbb. En este caso obtendremos un mapa de la Ciudad Autónoma de Buenos Aires, ya que vamos a localizar las Postas de Vacunación que geocodificamos en el paso anterior.

bbox <- getbb("Ciudad Autónoma de Buenos Aires,Argentina")

head(bbox)
##         min       max
## x -58.53145 -58.33514
## y -34.70564 -34.52655

Ya tenemos los límites de la ciudad en la variable bbox. Ahora vamos a descargar la capa base utilizando ggmap. Lo que permite hacer ggmap es, en esencia, es añadir a los gráficos ya conocidos una capa cartográfica adicional. Para eso usa recursos disponibles en la web a través de APIs; ggmap tiene funciones para obtener mapas (de diversos tipos y de distintos orígenes: Google, Stamen, OpenStreetMap). Nosotros vamos a obtener un mapa de Stamen que es un estudio de visualización de datos que abre sus mapas para el uso de la comunidad. Pueden ver sus diferentes diseños en Stamen maps

#Vamos a descargar el mapa

CABATerra <- get_stamenmap(bbox = bbox, #En bboxle indicamos los límites específicos, ya lo sabemos!
                      maptype = "terrain", #En maptype le indicamos el tipo de mapa que queremos descargarnos. 
                      zoom=13)#En zoom le indicamos el nivel de detalle que queremos. A más zoom, más detalle, pero también más pesado el archivo que nos descarguemos…

#que tenemos en la variable CABA?
head(CABATerra) #raro?
## [1] "#EAEAEA" "#FAFAFA" "#D1D1D1" "#999999" "#999999" "#999999"
CABATerra #R ya nos dice en la descripción que estamos trabajando con una imagen de mapa de Stamen Maps... e incluso nos sugiere como graficarlo!!!
## 1268x1143 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

Vamos a hacerle caso a los consejos de R y vamos a utilizar ggmap para dibujar el mapa. Por suerte para nosotros, ggmap tiene la misma sintaxis que ggplot(), por lo tanto no tendremos demasiados inconvenientes.

ggmap(CABATerra)

Y si probamos con otro tipo de mapa?

CABATonner <- get_stamenmap(bbox = bbox, 
                      maptype = "toner-background",  
                      zoom=13)

ggmap(CABATonner)

Excelente, tenemos nuestra capa de base, pero ahora cómo ubicamos en contexto las Postas de Vacunación de la Ciudad? Fácil… agregando una nueva capa a nuestro mapa!

ggmap(CABATonner) +
  geom_sf(data=Postas_Vacunacion_Direccion, inherit.aes=FALSE, color = "red", size = 2)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

¿Qué les parece el mapa que obtuvimos?…